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Abstract 


A numerical method is presented for solving the artificial compressibility form of the two-di- 
mensional time-dependent incompressible Euler equations. The approach is based on using an 
approximate Riemann solver for the cell face numerical flux of a finite volume discretization. Char- 
acteristic variable boundary conditions are developed and presented for all boundaries and in-flow 
out-flow situations. The system of algebraic equations is solved using the discretized Newton-re- 
laxation (DNR) implicit method. Numerical results are presented for both steady and unsteady flow. 
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I. Introduction 


As a predecessor to the three-dimensional time-dependent dynamic multiblock incompressible 
laminar and/or turbulent Navier-Stokes work [1 and 2], the two-dimensional incompressible Euler 
equations were investigated and a research code was written to solve the equations numerically. It 
turns out that this rather simple two-dimensional incompressible Euler code has been extremely use- 
ful in investigating various numerical flux formulations, boundary conditions, and numerical solu- 
tion algorithms that have fed back into the three-dimensional codes. It has also proven to be a useful 
teaching tool for introducing computational techniques to various interested parties. Moreover, it 
has been used as a test code for parallelization. In view of a number of requests for information, 
it seemed appropriate to record the details of the particular numerical solution scheme used to solve 
these time-dependent incompressible Euler equations. 

The computer code resulting from this effort is part of a family of Mississippi State codes that 
are now known as UNCLE (Unsteady Computation of Field Equations). Since this code is used 
as a research tool in solving the equations of Leonhard Euler, it is named UNCLE.LEONH ARD. 

The particular numerical methods used to solve the equations are the current ones used as of this 
writing. No attempt is made to review all of the techniques that have been tried and used, and/or 
tried and discarded. The emphasis is on simplicity, while simultaneously providing information on 
numerical methods that form much of the foundation of the three-dimensional time-dependent 
Navier-Stokes solvers. 

The system of equations to be solved is based on the artificial compressibility idea of Chorin[3]. 
This nonlinear system and its mathematically equivalent quasilinear system is presented in Chapter 
II. The eigensystem of the flux Jacobian matrices of the quasilinear system is developed in Chapter 
TIT A finite volume discretization is used and this formulation is given in Chapter IV. The cell face 
numerical flux is based on an approximate Riemann solver which uses the eigensystem of the quasi- 
linear form of the equations and this flux formulation is developed in Chapter V. The treatment of 
boundary conditions is given in Chapter VI. The numerical method used to solve the resulting sys- 
tem of algebraic equations is developed and presented in Chapter VII. Numerical results are pres- 
ented in Chapter VTII and concluding remarks are offered in Chapter IX. 

The body of this report is concerned with time-dependent equations on a stationary grid. How- 
ever, the class of problems that can be solved is expanded considerably if the grid is allowed to move. 
Therefore, since the basic methods presented can be applied to dynamic grids, an appendix is in- 
cluded that contains the two-dimensional time-dependent incompressible Euler equations on a dy- 
namic grid, along with the corresponding eigensystem. 
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II. Equations 


The two-dimensional time-dependent incompressible Euler equations in Cartesian coordinates 


are 


*1 + K + «£ = o 

dt dx dy 


( 2 . 1 ) 


where 


9 = 


/ = 


0 
U 
V 

u 

u 2 + p 

uv 


8 = 


uv 


v l + p 


These equations were made dimensionless by the relations 
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( 2 . 2 ) 


where the barred quantities are dimensional variables, L is a characteristic length, q r is the magni- 
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tude of a reference velocity, g r is a reference density, and t is time. The density is assumed to be 
constant and hence q is unity. The velocity components u and v correspond to the x and y directions 
respectively. Notice that the pressure variable p is twice the pressure coefficient. 

The three-dimensional time-dependent incompressible Euler equations in Cartesian coordi- 
nates were transformed to time-dependent curvilinear coordinates in Ref. 1 . After the equations 
were transformed to curvilinear coordinates they were written in Chorin’s [3] artificial compressibil- 
ity form as explained in Ref. 1 . Following this same approach, but for the more simple case of two 
dimensions rather than three and for a stationary curvilinear coordinate system as opposed to one 
that can move in time, the artificial compressibility form of the two-dimensional time-dependent 
incompressible Euler equations in curvilinear coordinates is 

dQ + dF , dG = ft (2.3) 

dr d£ drj 

where 
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lx = j~ x y v 

!y = - J~ l X V 

Vx = - J~ l 
T)y = J~ l X £ 

and fi = constant. Notice that flux vectors F and G can be written 


where 


K = J 


pe k 


ud k + k x p 


vd k + k y p 


8 k = k x u + k y v 


(2.4) 


(2.5) 


K = F and 6 k = U for k = £ 
K = G and 6 k = V for k = rj 


The quasilinear form of Eq. (2.3) is 


where 


d S + A§ + = o 

dr dr) 


A 


dF_ d dG_ 

9Q • se 


( 2 . 6 ) 
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The Jacobian matrices A and B are denoted by K where the elements of K are given by 


K« 


dK t 

dQj 


Therefore, 


and 


K — F and K = A when k = £ 


K = G and K = B when k = rj 


( 2 . 8 ) 


The flux vectors represented by Eq. (2.4) can be written in terms of the elements of the dependent 
variable vector Q = [ Jp, Ju, Jv] T = Q x Qi, Q3] as 


K = 


(k x Q2 + ky Q3 ) 
j 0-2 ky Q 3 )] F k x ( 2 j 

j Q 2 ky < 2 3 )] + k y Q\ 


(2.9) 


The elements of the Jacobian matrix K are determined from Eq. (2.8) resulting in 


K = 


0 
k r 
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( 2 . 10 ) 
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III. Eigensystem 


The development of the numerical method used to solve the Euler equations (Eq.(2.3)) depends 
on the eigenvalues and eigenvectors of the flux Jacobian matrix K given by Eq. (2. 10). These eigen- 
values and eigenvectors comprise the eigensystem of TC, which is developed in this Chapter. 

The eigenvalues of K are given by the solution to the characteristic polynomial resulting from 

I K - X I I =0 (3.1) 

It is relatively easy to solve this polynomial and obtain the eigenvalues 


<£> 

II 

(3.2a) 
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II 

<N 

(3.2b) 

I 

II 

(3.2c) 

where 


c = dj + p {kl + kj)J 

(3.3) 

and 6 k is given by Eq. (2.5). 


Because the eigenvalues given by Eqs. (3.2) are real and distinct, the system is strictly hyperbolic 
and thus one is assured of a linearly independent set of eigenvectors [4]. The right eigenvectors, r, 
are obtained from the solution to 

[K - X I) r = 0 

(3.4) 


A set of linearly independent right eigenvectors is given by 

0 
k y 

kx 
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A similarity matrix R k can be formed by using these right eigenvectors as columns 


Rk = 
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— k v 


uX 2 

T 

vX 7 
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+ k x 


+ k v 


uX 3 

7" + kx 

v h. + k 
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(3.6) 


A set of left eigenvectors, y T , can be obtained from the solution to 


y T (K - A 7) =0 


(3.7) 


In this case, since the R k matrix is only 3x3, it is probably easier to invert R k and obtain 
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where the rows of R k 1 are the left eigenvectors. 

The Jacobian matrices are now diagonalized by 

R~ l KR k = A k (3.9) 

where A k is the diagonal matrix with diagonal elements X x ,X 2 , and/t 3 . That is 
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0 


0 



where the subscript it corresponds to k — £ or k — rj, and 0 k is given by Eq.(2.5) and c is given 
by Eq.(3.3). 

The eigensystem of K is therefore composed of the eigenvalues given by Eqs.(3.2), the right ei- 
genvectors given as the columns of the matrix R k (Eq.(3.6)), and the left eigenvectors given as the 

rows of the matrix R k 1 (Eq.(3.8)). Note that each eigenvector corresponds to one particular eigen- 
value. 
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IV. Discretization 


MacCormack and Paullay [5] discretized the equations into what is now commonly referred to 
as the finite volume form. By integrating Eq. (2.3) in space with the index notation shown in Fig. 
1, the finite volume form of the equations can be written 


dQ 

dr 


mr < + 



(4.1) 


where the bar over the time derivative indicates various numerical approximations to the time rate 
of change of the dependent variable vector, which is recognized in this approach to be located at the 
center of the cell. Using central difference notation, Eq. (4. 1) can be written 


dQ 

dr 
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AS 
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VG 

Arj 


= 0 


(4.2) 


With A£ = 1 and Arj = 1, Eq. (4.2) becomes 


dQ 

dr 


) + <5, F + d. G = 0 

i,j 


(4.3) 


Various approximations of the time derivative have been used. The most commonly used approxi- 
mation is first-order in time given by 


dQ 

dr 
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(4.4) 
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where the superscripts indicate time level. For higher order in time the second-order backward- 
approximation 



71—1 


3Q n + l -4 Q n . + Q n 
Zhl Zhl 
2 Ax 


(4.5) 


is commonly used. 

The finite volume formulation given by Eq. (4.3) states the relation between the summation of 
fluxes across the faces of a cell and the time rate of change of the dependent variables in the cell. 
The finite volume formulation offers considerable flexibility with regard to geometric complexity 
because the cells can be of any shape in physical space. 


Use of central difference notation in this finite volume formulation (Eq. (4.3)) simply indicates 
the difference in flux vectors of opposite faces for the cell centered at the point i, j and should not 
be construed to mean the numerical solution scheme is central difference. The reason the scheme 
will not be central difference is that the flux vectors will be determined using upwind (one-sided) 
or upwind— biased techniques as explained in the following Chapter. 


10 


V. Numerical Flux 


In this finite volume discretization the summation of fluxes passing across cell faces is balanced 
with the time rate of change of the dependent variables in the cell. There are numerous ways of ob 
taining the flux vector at cell faces; however, the quality of the numerical solution is critically depen- 
dent on how the flux vector is formulated. Having investigated numerous ways of forming the flux, 
it is the authors’ observation that Godunov[6] type schemes have much to offer. In this approach 
a series of Riemann problems are solved where each Riemann problem corresponds to a cell face. 
Roe [7] followed Godunov’s idea, but developed a much more computationally efficient method of 
solving the Riemann problems. This class of numerical approaches is referred to as approximate 
Riemann solvers and are frequently used in the computation of compressible flow where the convec- 
tive portion of the equations form a hyperbolic system. By using Chorin’s [3] idea of artificial com- 
pressibility, the incompressible equations are also hyperbolic and therefore approximate Riemann 
solvers can be used. The approximate Riemann solver of Roe is the foundation of the present numer- 
ical approach. 

A tutorial type discussion of the Roe scheme is given in Ref. [8]. In the development of this 
method, waves are assumed to move normal to a cell face. This allows the numerical flux vector 
to be obtained independently in each of the computational coordinate directions £ and tj . It is popu- 
lar in the literature to state three different ways of getting the flux vector at a cell face. The most 
frequently used method, usually written as the third expression for the flux vector, is actually an aver- 
age of the first two expressions [7]. There is a larger operation count in using this averaged expres- 
sion as compared to either one of the other two, and so the expression used here for the first-order 

Roe numerical flux vector at cell face i + in the £ computational direction is 



(5.1) 


where 


a®, .= <®, . 

i 2 yj * ■ *^* 2 ’ J 


[Qi+\,j Qi ,j) 


(5.2) 
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and X correspond to the negative eigenvalues of the Roe matrix, r® are the right eigenvectors 
of the Roe matrix, P are left eigenvectors of the Roe matrix, and the scalars are jumps in the 
characteristic variables. The subscript i + ±,j in the equations above indicates that the metrics used 

correspond to the cell face located at / + The dependent variables in the eigenvalues and eigen- 

vectors are the Roe averaged variables [7] at cell face i + ^ j since they comprise the eigensystem 
of the Roe matrix. The flux vector F(Q t J) in Eq. (5.1) is evaluated using. the dependent variable 

vector as indicated (not the Roe averaged variables), but the subscript / + ^,j on the bracket indi- 
cates the metrics at / + ^ ,j are to be used. The eigensystem used in Eqs. (5. 1) and (5.2) are, there- 
fore, the columns of Eq. (3.6) for the right eigenvectors, the rows of Eq. (3.8) for the left eigenvec- 
tors, and the diagonal elements of Eq. (3.10) (with the notation that Xj in Eq. (3.10) is X® in Eq. 
(5.1)), all with the Roe averaged velocities 

u i+\,j = 2 ( M ' + 1 ’•/ + “'••/■) (5.3a) 

and 

Vi.; = £(v 1+ i, , + (5.3b) 

and with it = £ in Eqs. (3.6), (3.8), and (3.10). The flux vector F{Q t j) in Eq. (5.1) comes from 
Eq. (2.4) with k = £ and is 

» Hj 

+ Pi’i 
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where 



U iJ = u i.j + ^ v /j 


(5.5) 


Another way of writing Eq. (5.1) is [8] 




+ A '(Q,.j • e i+ ,j)(fi (+ ,j 



(5.6) 


where A is the Roe matrix given by 


A = R^ 1 (5-7) 

where A| has only nonpositive elements. All the dependent variables in Eq. (5.7) are the Roe aver- 
aged variables given by Eqs. (5.3). All the metrics in the A matrix correspond to cell face / f + i ,j. 

The numerical flux vector, G, for the rj direction can be written in a fashion analogous to the 
vector F for the £ direction. For example, the vector G in the form similar to Eq. (5.6) would be 


e, J+i = [ c(a,j) 


' irj+i 


+ B (Q itJ , Q iJ+l )(Q iJ+l ~ Q itJ ) (5.8) 


where B is the Roe matrix given by 


B — Rfj Ajj Ryj 


(5.9) 
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where A ^ has only nonpositive elements. The dependent variables in Eq. (5.9) are the Roe aver- 
aged variables given in the r\ direction by 


*<j+r2 ( v m+i + v w) 


(5.10a) 


(5.10b) 


All the metrics in the B matrix correspond to cell face i ,j + i. 

Consider the numerical flux vector in the £ direction given by Eq. (5.6). Notice that since Eq. 
(5.6) is the solution to a Riemann problem, the conditions to the right of the cell face at i + i ,j are 
given by Q i+ , ; and the conditions to the left of this cell face are given by £>, j. Equation (5.6), 
therefore, could be written as 




+ A 



(5.11) 


where 


er + i,-e i+1J 


and 


- a j 


The corresponding Roe averaged variables used in A would be 


«. + 1 .-if «* .+ <, •) 


(5.12a) 


(5.12b) 


(5.13a) 
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and 


7 iHj = j{ v U 


(5.13b) 


The numerical flux vectors introduced thus far lead to first— order accurate schemes in space. 
These numerical schemes can be made higher order by using the MUSCL approach of van Leer [9]. 
Anderson, Thomas, and van Leer [10] used a MUSCL-type approach in flux vector split schemes. 

ft 

Following [10], the dependent variables just to the right of the cell face, Q , located at 
i + 1 J and just to the left of the same cell face, Q L , are written as 


= G. + I .) - * [ (1 - *)(&«,,• - C<+i,;) + U + x)(2. + i .i - &.;) ] (5.14a) 


= Qi.t + | [ l 1 - *iQtj ~ Qi-u) + (* + *)(a + i j - Qi.i) ] 


(5.14b) 


For 0 = 0, Eqs. (5.14) recover Eqs. (5.12) and hence the numerical scheme would be first-order 
in space. For higher order schemes set 0 = 1 . With x = - 1 only points to the right of the cell 

face are used for Q R and only points to the left of the cell face are used for Q , and a second-order 
scheme results. With x = 1/3 two points to the right and one point to the left of the cell face are 

used for Q R , and two points to the left and one point to the right of the cell face are used for Q L . 
The resulting scheme is referred to as third-order upwind-biased [10]. Note that the third-order 
upwind-biased scheme depends on information from the same number of points as the second-order 
scheme. 

The same approach as used here for the F numerical flux vector corresponding to the £, or i, 
computational direction is applied to the G numerical flux vector corresponding to the rj, or j, com- 
putational direction. For example, G would be 
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Note from Eqs. (3.2) and (3.3) that for incompressible flow, unlike compressible flow, the se- 
cond eigenvalue is never less than zero. Note, also, that the numerical flux vector formulation used 
for F and G in Eqs. (5.6) and (5.8) do not use the positive eigenvalues. Advantage is taken of this 
fact during coding in order to reduce memory and floating point operations. To illustrate this, con- 
sider, for example, the second term on the right-hand-side of Eq. (5.6) ( or Eq. (5.1) since they are 
the same) and use the following generic definitions 


A = R A R 


(5.16) 


where 


'll 

r \2 

r 13 
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r 33 
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Then 


u 


A • dQ = Ajf/jj 6Q X + l l2 &0-2 + l \3 ^Q-i) 


'21 


31 


12 


+ A 2 (/ 2 1 + l 22 &Qi + ^23 ^ 3 ) 


r 22 


'32 


rni 


+ A 3 (/ 31 6Q x + l 32 6Q 2 + Z 33 dQ 2 ) 


r 23 

r 33 


(5.17) 


Because X 2 will never be less than zero the second term in Eq. (5.17) will never be needed, which 
means that the second left and right eigenvectors will never be needed. Moreover, the first term will 
be needed only part of the time, i.e. when X x < 0. 

It should be pointed out that the alternative flux vector formulation of Roe [7] given by 




e i+ i j)(e i+ , j - a.,) 


(5.18) 


where 


_+ 

A 


= R t A 1 
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and Ag has only non-negative elements, could have been used. Because X 3 will never be greater 
than zero (see Eq. (3.2c)), either formulation would lead to the same savings and, more importantly, 
to the same numerical results. 

It might be remarked that the approach just outlined for obtaining the numerical flux vector is 
a significant departure from the approach used in Refs. [1] and [2]. The reason is that the method 
previously used came from experiences gained in compressible flows [8] and [11]. Limiters were 
used in compressible flow and the method used to formulate the numerical flux vector in [8] was 
found to have certain advantages compared to a MUSCL-type approach. However, limiters have 
never been required in incompressible flow and the approach just outlined, and the approach of [8], 
tend to give the same numerical result for incompressible flow as will be demonstrated in Chapter 
VTII. Since the convergence rate of the two methods has been found to be essentially the same, and 
because the present approach requires less memory and less floating point operations, the present 
approach is currently the preferred approach. 
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VI. Boundary Conditions 


Boundary conditions on all surfaces are developed using characteristic variable boundary condi- 
tions. The approach follows that presented in Ref. [12] for the numerical solution of the three-di- 
mensional unsteady Euler equations for compressible flow. Loper [13] refers to this approach as 
“reference plane” boundary conditions, where the reference plane corresponds to the plane resulting 
from keeping constant the computational coordinate along which the boundary conditions are devel- 
oped. In what follows, this boundary plane will be located on the plane resulting from keeping k 
= constant, where k can be either the £ or rj computational coordinate. Then from Ref. [12], the 
characteristic variable boundary conditions will depend on the equation 


l (d_Q 

k \dT 


+ R k A k R 


-i dQ\ 
k dk) 


= 0 


( 6 . 1 ) 


or 




( 6 . 2 ) 


Assuming R k 1 to be a constant matrix with elements corresponding to some reference conditions 
denoted by the subscript “o”, Eq. (6.2) can be written 


dW, 

dr 


+ A, 


Wk 

dk 


= 0 


(6.3) 


where W k is the characteristic variable vector and is defined as 

W k = R~' a Q (6.4) 

Because A k is a diagonal matrix, Eq. (6.3) represents three uncoupled equations. For example, 
if the elements of W k are given by 



then Eq. (6.3) can be written as 


dw k,m 

dr 


+ k m 


d w k,m 

dk 


= 0 


( 6 . 6 ) 


for m=l,2, and 3. Note that Eq. (6.6) can also be written as 


dw k,m _ dw k,m , dW k,_m_ n (6.7) 

dr dr + Am dk 


where 


X 


m 


= dk 

dr 


( 6 . 8 ) 


Therefore, each characteristic variable is associated with one particular eigenvalue, and this charac- 
teristic variable is constant along the characteristic line described by the direction given by Eq. (6.8). 

Since the left eigenvectors are rows of the matrix R k 1 , each characteristic variable w k _ m is the 
dot product of the m‘ h left eigenvector with the dependent variable vector Q. Therefore, the charac- 
teristic variables are obtained using R k 1 given by Eq. (3.8) evaluated at some reference conditions 
times Q. However, note from Eq. (6.7) that any constant times a characteristic variable also satisfies 
Eq. (6.7). This fact can be used to simplify the characteristic variables somewhat. Taking advantage 
of this fact is helpful in order to simplify the algebra involved in developing explicit expressions for 
the boundary conditions. The characteristic variables used here are 


w k t = (u k y — v k x ) 0 p — (v 6 k + (3 ky) 0 u + (w 0 k + k x ) 0 v (6.9a) 

w kt2 = {6 k - c)o p - P(k x u + ky v) (6.9b) 

w k 2 = (fi/c "b c)<? P $(.k x u + k y v) (6.9c) 


where the characteristic variables correspond in order to the eigenvalues 
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A 2 

*3 

As a reminder 


= 0* 

(6.10a) 

II 

< 3 ^ 

+ 

O 

(6.10b) 

1 

<£> 

II 

(6.10c) 


and 


6 k = k x U + ky V 


( 6 . 11 ) 



+ + 



( 6 . 12 ) 


These equations are now used to develop far field and impermeable wall boundary conditions for 
the two-dimensional incompressible Euler equations. 


Codirectional Inflow: 

In the development of boundary conditions for the compressible Euler equations on dynamic 
grids, Janus [14] coined the terminology codirectional and contradirectional. This can be a conve- 
nient terminology, and it is used here for the incompressible Euler equations on stationary grids. 
From Fig. 2a, codirectional flow is shown to be flow traveling across the boundary of interest in the 
direction of increasing k computational coordinate. Therefore, flow coming into the computational 
domain traveling in the direction of increasing k is codirectional inflow. In Fig. 2 the point labeled 
a corresponds to flow approaching the boundary, the point “b ” corresponds to the boundary, and 
the point labeled corresponds to flow leaving the boundary. 

For codirectional inflow the first two eigenvalues of Eqs. (6.10) are positive and run from out- 
side the computational domain toward the boundary. The third eigenvalue is negative and runs from 
inside the computational domain toward the boundary. Since the respective characteristic variables 
along these characteristics are constant, reference to Fig. 2a indicates that 
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5 

Cs 

II 

"i" 

(6.13a) 

M a = K 4 

(6.13b) 

W = 

(6.13c) 


That is, the first two characteristic variables are constant from point “a ” to point “b ” and the third 
characteristic variable is constant from point “/’’to point “b”. Point “a ” is taken as some reference 
point. For example, for flow about an airfoil in free air the characteristic variable is evaluated using 
freestream conditions. The point “l” is taken as the center of the first cell inside the computational 
domain. This leads to the equations 


(uky-V k x ) o Pb - (v e k + p ky) o U b + (u e k + p k x ) o V b = 
[u ky — v k x ) o — (v 6 k + ft k y ) o u oo + (w 6 k + (i k x ) o v«, 


(6.14a) 


( 6 k - c ) a Pb - P k x u b - P k y v b = ( e k - c ) 0 Pa* - P °k,c o (6.14b) 

( d k + c ) 0 Pb ~ P k x u b ~ pk y v b = (0 k + c) o Pl -p O k l 


Equations (6.14) can be solved for the three unknowns p b ,u b , and v b to obtain 


Poo + P/ 
2 

Ub ~ U co 
V b = V oo 


2~ ®k,o (p°° Pi) ] 

(6.15a) 

p Co b “ C °) + ^ kx ] 

(6.15b) 

Pc 0 b [ V °( 6 ^o ~ c o) +P k y] 

(6.15c) 


The reference conditions indicated by the subscript “o ” are constructed from an average of de- 
pendent variables located at the boundary and at the center of the first cell inside the computational 
domain at the previous time step. 
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Contradirectional Inflow: 

Contradirectional inflow is illustrated in Fig. 2b. In this case the first and third eigenvalues given 
by Eqs. (6.10) are negative and run from outside the computational domain to the boundary. The 
second eigenvalue is positive and runs from inside the computational domain to the boundary. The 


characteristic variables are, therefore, grouped as 



KO 

Cs 

II 

'IT 

*- 


(6.16a) 

K 2, 

1, = K4 


(6.16b) 

[ w k, 3 ) 

L = W ( 


(6.16c) 


These relations lead to the equations 


(uk y -v k^ o p b -(v d k + p ky) o u b + (u9 k + p k x ) o v b = 

ky V o p OO ^ V 0 k kyj o M oo ""t" 9 k T kjcj V oo 


(6.17a) 


( e k - c ) 0 Pb - P k x u b - P k y v b = ( 6 k - c ) 0 Pi ~P k * «/ - P k y V/ (6.17b) 


(' 9 k + c) o p b - (ik x u b - ftky v b = ( 9 k + c) o p oo - fik x u oo - ($k y Voo (6.17c) 


Equations (6.17) can be solved for p b ,u b , and v b to obtain 


Pb = 


Pec + Pl 
2 



@k , /) &k,o ip a Pi) ] 


(6.18a) 


U b = Moo 



(6.18b) 


v b = 


P °° P b 

fic 0 




(6.18c) 


23 


Codirectional Outflow: 

In this case the first two eigenvalues are positive and run from inside the computational domain 
to the boundary. The third eigenvalue is negative and runs from outside the computational domain 
to the boundary. However, unlike inflow, the characteristic variable associated with this third eigen- 
value is not easy (or perhaps not enough is known about it) to evaluate. Therefore, for outflow, the 
pressure is specified on the boundary, and this specification of pressure eliminates, or replaces, the 
relation involving the incoming characteristic variable. The equations to solve for u b and v b , there- 
fore, are 


( v e k + P k y) 0 u b + ( M 6 k + P k x) a V b = (uk y -v k x ) o (p a - p b ) 

- { v e k + P k y) 0 u a + ( M e k + P k x) a v a (6.19a) 

- P k x Ub - P k y V b = {&k ~ c ) 0 ( Pa - Pb) - P k x Ua - P k y V a (6.19b) 


Pb P exit 


(6.19c) 


Recall that the point “a ” (see Fig. 2a) is associated with approaching flow and in this case it is located 
inside the computational domain at the center of the first cell inside the domain. With p b specified 
by Eq. (6.19c), Eqs. (6.19a) and (6.19b) can be solved for u b and v b yielding 


Pb = Pexit 


(6.20a) 


[ “oftlc.o ~ c o) +P k x 

(6.20b) 

V ° + Pc 0 

[ V o[@k,o ~ c o) + P k y] 

(6.20c) 


Contradirectional Outflow: 

In this case the first and third eigenvalues are negative and run from inside the computational 
domain to the boundary. The second eigenvalue is positive and runs from outside the computational 
domain to the boundary. The characteristic variable relation corresponding to this incoming charac- 
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teristic is replaced by the specification of pressure on the boundary just as was done for the codirec- 
tional outflow case. The equations for the boundary conditions, therefore, are 


( v 6 k + P k y) 0 u b + ( u Ok + P k x) 0 v* = (u k y - v k x ) o (p a - p b ) 

- (v d k + k y ) o u a + (u 6 k + ft k x ) o v a (6.21a) 


Pb ~ P exit 


(6.21b) 


- pk x u b - Pky v b = (e k + c) o (p a - p b ) - Pk x u a - pky v a (6.21c) 


These equations can be solved to give 


Pb Pexit 

(6.22a) 

pc 0 [ U °( e k,o + Co) + P k x j 

(6.22b) 

p Co [ V °{®k , o + c o ) + P *>•] 

(6.22c) 


Impermeable Surface: 

Consider an impermeable surface corresponding to a constant k line where the computational 
domain is located in the direction of increasing k. Then the center of the first cell inside the computa- 
tional domain can be connected to the impermeable surface by a negative running characteristic cor- 
responding to the negative eigenvalue. The characteristic variable relation corresponding to the pos- 
itive eigenvalue is then replaced by the physical fact that there is no flow through the impermeable 
surface, or 


Ot,b = 0 


( 6 . 23 ) 


If the computational domain is in the direction of decreasing k from the impermeable surface, the 
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center of the first cell inside the computational domain can be connected to the impermeable surface 
by a positive running characteristic and Eq. (6.23) replaces the characteristic variable relation corre- 
sponding the negative eigenvalue. This leads to the following characteristic variable relationship 


Co P b *p Qk,b = C o Pin T P Qk,in 


(6.24) 


The 6 k 0 term in Eqs. (6.9) has been set to zero in arriving at Eq. (6.24), and the subscript "in” in 
Eq. (6.24) identifies the center of the first cell inside the computational domain. By using Eq. (6.23) 
in Eq. (6.24) an expression for the pressure at an impermeable surface can be obtained as 


Pb Pin 



(6.25) 


where the minus sign corresponds to the situation where the computational domain is in the positive 
k direction relative to the boundary, and the plus sign corresponds to the computational domain being 
in the negative k direction relative to the boundary. 

Pressure is formally the only thing that is needed at an impermeable boundary for an Euler equa- 
tion solver. However, to use a flux formulation higher than first order at the first cell face from a 
boundary, the remaining dependent variables are needed at the surface. An approximate expression 
for these remaining variables can be obtained by using the first characteristic variable (Eq. (6.9a)). 
The first eigenvalue is zero at an impermeable surface and consequently the first eigenvalue does 
not point from inside the computational domain to the surface. However, if the approximation is 
made that the first characteristic variable relation can be written as 


[u ky- V k x ) o p b - pky u b + fik x v b = (u ky- v k x ) o p ^ - fik y u in + fik x v in 


(6.26) 


then both additional velocity components can be obtained by solving Eqs. (6.26) and (6.23) for 
u b and v b . The complete impermeable surface boundary conditions can be written as 
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Pb - Pin T 


(6.27a) 


P_ @k,in 
c o 


u b = u in + 


Va = v,„ + 




k x | 

C ° k l + k y) 

Vo _ k y \ 

C ° k l + k y) 


®k,in 

®k,in 


(6.27b) 

(6.27c) 


where, again, it is noted that the minus sign corresponds to the computational domain being located 
in the positive k direction relative to the boundary, and the plus sign corresponds to the computational 
domain being located in the negative k direction relative to the boundary. 
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VII. Numerical Solution 


7.1 Formulation 


An implicit solution scheme is developed by using the numerical flux vectors F and G from 
Eqs. (5.11) and (5.15) in the finite volume descretization given by Eq. (4.3) to obtain 

AO n 

-Zhi + d,F(Q"*') + d j G(Q" + ') = 0 (7-1) 

where 



AQij = Qij ~ 

eh 

(7.2a) 

<5, F(0 n+1 ) 

= J W e " +l ) 

-r.-i.dO 

(7.2b) 

djG(Q n+l ) 

-wo 

-wo 

(7.2c) 


The first-order time derivative given by Eq. (4.4) is used in Eq. (7. 1) for illustration. The approach 
used to develop methods of higher-order accuracy in time is the same as that presented below for 
first-order time accuracy. 

The normal procedure for the solution of Eq. (7.1) would be to linearize the spatial difference 
terms, move the terms not containing A Q n to the right-hand-side of the equations, and solve for 
A (2". This is particularly true for problems expected to have steady state solutions, because the sum 
of the spatial difference operator terms as well as A Q n would both go to zero. However, for unsteady 

flow, the ideal situation would be to find Q n + 1 such that Eq. (7.1) is satisfied. That is, find Q n + 1 
such that [15] 


T(e” +1 ) = o 


(7.3) 


where 
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(7.4) 


v(Q n+l ) = ^r + d i F (Q n+X ) +djG(Q n+ ') 

One way of attempting to solve this problem is to use Newton’s method. Newton’s method [16] 
for the function £F(x) would be 


3 : '(x m )(x m+1 -x m )= - 


(7.5) 


where m = 1,2, 3, ... and T'(x) is the Jacobian matrix of the vector T(x). With Q n+ 1 = x, Eq. (7.5) 
becomes 

+ 1, m + 1 r \n + 1, m 


9 r '(Q n+1 ’ m )(<2" 


) = - 9 : (0 n+1 ’ m ) (7.6) 


In principle the generated sequence Q n + 1 ’ m + 1 converges to Q n + 1 and, hence, Eq. (7.3) is satisfied. 

Newton’s method has the attractive feature that the convergence rate can be quadratic [16]. 
Whereas, quadratic convergence is possible for computational fluid dynamic problems it is difficult 
to achieve; and even if achieved, it may not be the most computationally efficient way to solve the 
problem [17]. It usually occurs in practice that in the interest of computer time and memory one 
does not strive for quadratic convergence, rather one seeks variations [16] of Newton’s method. Be- 
cause the numerical problem under investigation is well formulated in the context of Newton’s meth- 
od, particularly for unsteady flows since the residual (which is driven to zero ) now contains the time 
derivative, this approach of using variations or approximations of Newton’s method is followed 
here. 


7.2 Jacobian 

The solution to Eq. (7.6) formally requires the Jacobian matrix 9F'(x), where the elements of 
T'(x) are defined as 


a ij (*) = 


d^i (x) 


dxj 


(7.7) 
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Analytical development of the Jacobian of the Roe numerical flux vector is not simple, and evidently 
impractical to obtain in three dimensions. Since it is the intent of this work to produce techniques 
that can be used in three dimensions, an analytical derivation of the Jacobian doesn’t appear feasible. 
In some cases flux vector splitting [18] has been used to obtain approximate analytical Jacobians 
for the solution matrix while flux differencing splitting, such as Roe’s method, has been used on the 
right-hand-side of the equations [19]. Whereas this approach is inconsistent in that flux vector split- 
ting is used on the left, and flux difference splitting is used on the right, it has been shown to work 
quite well [19 and 20]. This approach, however, was used in compressible flow where the flux vector 
was homogeneous of degree one and splitting of the flux vector was straightforward. For this incom- 
pressible formulation the flux vector does not have this homogeneity property and it is not obvious 
how to follow this approach. Therefore, the approach used here is to obtain the Jacobian numerically 
by using difference quotients [11]. The Jacobian elements are obtained from 


a i j (*) = 


SF,- (x + h ej) — T, (jc) 
h 


(7.8) 


where ej is the j th unit vector and 


h ~ J machine e 


(7.9) 


When a numerical Jacobian is used, Eq. (7.6) is known as discretized Newton iteration [16]. Define 
the following 


DF 


dF . , . 
1 + 2 './ 


i + 2’j dQi j 


). 


dG. .. i 


DG , .= ^ 

J + 2’J dQi j 'm 


‘■ >n ). 


(7.10a) 

(7.10b) 


where the first subscript of DF and DG refers to the location of the cell face and the second sub- 
script refers to the location of the dependent variable vector that the numerical flux vector is differen- 
tiated with respect to. Using these definitions, Eq. (7.6) can be written 
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-BC . AQ"*':"'-DF : , . AQ n *''" 

J 2’J 1 l 'J 1 1,7 


+ 1, m 


(x 

Vdr 


+ (X + JDF 1 + DG , — DF. , . - ZJG. , .) J0" + 

(+ 2-' y+2-^ *~i.» y-j,y/ 


+ ^. i+1 ^^ I : ; + D5 )+iJ+1 z. e -;; 


^n+ 1, m 


g« + l,m _ 

A: - + <5, F(e" + '•“) + 5j g(q" + '-”) 


(7.11) 


where 


AQ n+hm = g" + l ,m+ i - Q n + hm (7.12) 

The term 7 fl is an identity matrix except the first diagonal element is zero in order to satisfy the true 
continuity equation. 

Note that although the higher order flux vectors given by Eqs. (5.11) and (5.15) involve the de- 
pendent variable vector at more than one point on either side of a cell face only one point on either 
side was used in Eq. (7. 11). This means that the Jacobian of the first-order flux vectors is used on 
the left, whereas, higher-order flux vectors are used on the right. There are two reasons for this. 
Firstly, by using more than one point on either side the band width of the solution matrix (the Jaco- 
bian) would increase, requiring more memory and more CPU time to solve the problem. Secondly, 
experiences gained in using more than one point have indicated that convergence is not improved 
a great deal when more points are used. For quadratic convergence, all points involved in the flux 
vectors should be included in the Jacobian. However, quadratic convergence is not obtained here 
and the approximation made by using only one point on either side of a cell face works rather well. 


During the course of this work, another formulation for the Jacobian in Eq. (7.11) was found to 
work rather well. In fact, it will be demonstrated that at least for the test case considered, it works 
better than Eq. (7.11). Note from the expressions for the higher-order numerical fluxes given by 

Eqs. (5. 11) and (5. 15), that they depend on Q R and Q L given, for example, by Eq. (5. 14). By form- 
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R Li R L 

ing Q and Q using Eq. (5. 14) and incrementing Q and Q by (h in the numerical derivative 

(Eq. (7.8)) one obtains an approximate numerical Jacobian that is computationally more efficient 
than using a strictly first-order numerical Jacobian. (That is, it is computationally more efficient 
when a higher-order flux is used on the right-hand-side of the equation.) The reason it is computa- 
tionally more efficient is that if a higher-order flux is used on the right in the residual term then this 
flux which has already been computed and is available, can be used in Eq. (7.8) and another first-or- 
der flux does not have to be computed as is necessary if a strictly first-order Jacobian is used on the 
left-hand-side. By using this slightly different numerical Jacobian, the system of equations to be 
solved that is analogous to Eq. (7.1 1) is given by 


a + 1 1 f n 7 ^-pr A + 1 , m 

DG. , _ AQ . . - DF. , r AQ t . 

J-\,L « tj— 1 i-\,L 


(-L 

\zlr 


+ DF i+lL + DG j+lL - DF i-lR 


DG H R ) AQ1+'-' 


+ DF 


i+i,R 


AQ 


n + 1 , m 

i+U 


+ DG 


j+\,R 


AQ 


n + l,m 

ij + 1 


V 


Q 


,n+ 1 ,m 
Fj 


~Q; ; 


— + <5, F(e" + , ’ m ) + dj G(Q n + 1 ’ m ) 


At 


(7.13) 


where the first subscript on DF and DG again refers to the location of the cell face and the second 
subscript refers to whether the numerical partial derivative is with respect to Q constructed from the 
right (R) or left (L). 


To begin the calculation at each time step, the following is used for the initial (m = 1) Q 


n + 1 , m 


Q" + U = G” 


(7.14) 


Note that when Eq. (7. 13) is solved for m = 1 only, the time derivative in the residual term is zero. 
This would correspond to the situation of linearizing only the spatial difference terms in Eq. (7.1). 
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This approach is used for steady state calculations and for some unsteady calculations when the vari- 
ations in time are not sever. This results in an approximate first-order time derivative solution. 

It was stated in Section 7.1 that the development of the solution process for higher-order in time 
is the same as for first-order in time. If the second-order backward time derivative given by Eq. 
(4.5) was used in place of the first-order time derivative given by Eq. (4.4), then the resulting system 
of equations analogous to Eq. (7.13) would be 


-DC , Jg n+I '" - DT. , ,AQ n *l' m 

j-x.L 1 ,-\,L .1 


+ [m + df ^.l + DG m.<. ~ df h,r - dg h,r) 


aq::; : '; + dg j+ , r aq^ 


3 Q n+hm - 40" . + Q n ~ 


n-\ 




7At 


(7.15) 


where, again, I a is an identity matrix except the first diagonal element is zero in order to satisfy the 

true continuity equation. Equation (7.14) is used in Eq. (7.15) for the initial value of Q n+] ’\just 
as was done for Eq. (7.13). Although it is not quite so easily seen, a solution to Eq. (7.15) for m = 
1 only, results in an approximate second-order time derivative solution analogous to the first-order 
time derivative solution resulting from Eq. (7. 1 3) when Eq. (7. 14) is used for an m = 1 only solution. 


Equation (7.13) for first-order in time or Eq. (7.15) for second-order in time constitute the sys- 
tem of equations that must be solved for the dependent variable vector Q n+ \ The solution of this 
system is the subject of the next Section. 


7.3 Solution of the System of Equations 

The numerical solution of Eq. (7.13) or (7.15) from time level n to time level n+1 calls for m 
iterations. The solution of these equations for each value of m is referred to as a Newton iteration. 
For unsteady flows, values of m from 1 to 5 are typically used. The case of m = 1 is sometimes 
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referred to as no Newton iterations since the equations solved are the same as those that would be 
solved if only the spatial difference terms of Eq. (7.1) were linearized. 

Note that the numerical solution of Eq. (7.13) or (7.15) for each Newton iteration amounts to 
solving a linear system of equations that is represented here as the linear algebra problem 

A x = b (7-16) 


i 

where x in this case is the change in Q from iteration m to m+1, A is the Jacobian at iteration 
m, and b is the residual vector at iteration m. Matrix A is a blocked banded matrix. The blocks are 
square matrices that are 3x3 in this case since there are three dependent variables at each point. Note 
from Eq. (7. 1 3) or (7. 1 5) that the Jacobian is evaluated at each m iteration. However, because x goes 
to zero (in principle) as m increases, the Jacobian has no influence on the converged solution so long 
as the Jacobian that is used will work. Therefore, in practice, the Jacobian is not updated at each 
Newton iteration in order to save CPU time. The solution process discussed below is appropriate 
whether or not the Jacobian A is updated after each Newton iteration. 


Write the matrix A in the form 


A = L + B + U 


(7.17) 


where L is a lower block triangular matrix with zeros on the diagonal, B is a block diagonal matrix, 
and U is an upper block triangular matrix with zeros on the diagonal. The first two terms of Eq. 
(7.13) or (7.15) contribute the elements of L, the third term contributes the elements of B, and the 
fourth and fifth terms contribute the elements of U. The equation to solve is 

(L + B+ U)x = b (7.18) 


Equation (7.18) is frequently solved in the CFD community by various methods of factoring, 
where the most commonly used method is probably that of Briley and McDonald [21]. The approach 
used here will be relaxation [22] . In many respects various factored schemes are similar to relaxation 
schemes, and some of these similarities are discussed in [11, 15, and 23]. 
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Equation (7. 18) is solved using the symmetric block Gauss-Seidel iterative method [24], This 
method performs a forward sweep through the computational domain by 


(L + B) * (1) + U x w = b 


( 0 ) 


(7.19) 


followed by a backward sweep 


L * (1) + (B + U) x {2) = b 


(7.20) 


This forward and backward iterative process can be written as 

(L + B) x (2p-1) + U x (2p ~ 2) = b 
L x (2p ~ l) + (B + U) x (2p) = b 

where p= 1,2, 3, ... Therefore, there are 2 p sweeps through the computational domain, half of 
them forward and half of them backward, for p complete symmetric block Gauss-Seidel iterations. 


(7.21a) 

(7.21b) 


To illustrate rather specifically the solution process, consider Eq. (7.21a) in the expanded form 


j x {2p ~ l) + L x {2p ~ l) + B- x {2p ~ l) 

L i J- 1 X iJ-l + L i~ 1 J X i- 1 ,j + >J i J 


4- II x^ 2p 2 ^ + U x^ 2p 2 ^ — b ■ 
+ U i+l J X i+\ j +U iJ+\ X i,j+\ D 1'J 


(7.22) 


In Eq. (7.22) each subscripted matrix L and U is a block ( 3x3 matrix in this case) and each subscripted 
vector x and b is a subvector (3 element vector in this case). Note that when one moves forward 
through the computational domain, then at an i,j location every subvector x behind this point has 
already been solved for and Eq. (7.22) can be written as 


35 



(7.23a) 


(2p— 1) r (2p-l) _ J (2p 1) 

B i,j x ij ~ j ~ Li j~\ X i j_ j i-l, y i-1, j 

(2p - 2 ) _ TT (2p-2) 

Ui+\,j X i+l,j ij+l i,j + 1 


In an analogous fashion, for the backward sweep given by Eq. (7.21b) one can solve the following 

, ( 2 P) 

equation for x 


(2 p) , j (2p — 1) t ( 2 P 0 

B ij - Kt ~ L ‘.i- ' x i j-i _ l <-' j V-M 


rj J2P) _ rj . , 

^«+l x i+l ,j l ’J +l ij + 1 


(7.23b) 


where p = 1, 2, 3, ... 

Although the overall solution process is an iterative method, the numerical solution of x at point 
i'j in Eqs. (7.23a) and (7.23b) is carried out using the direct method of Doolittle, which is a compact 
scheme for Gaussian elimination [25], Therefore, the overall process is iterative but the solution 
of the subsystems at each point ij is direct. 

In the numerical solution of the linear algebra problem given by Eq. (7. 16) it usually is impracti- 
cal to consider the solution in the form 


where the coefficient matrix A is inverted. The same is true here with regard to solving Eq. (7.16) 
because of the size of A involved in this problem. However, when a matrix is much smaller, like 
B is in the direct solution of Eq. (7.23a) and (7.23b), it may be of benefit to consider inverting the 
B block coefficient matrices. One is motivated to investigate this approach even more when it is 
noted that the Jacobian matrix, and hence the B block matrices, do not have to be updated for each 
Newton iteration, and certainly the B’s are not updated during the symmetric block Gauss-Seidel 
iterations. It turns out for the solution of these two-dimensional equations, it pays to invert the block 
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matrices B, which are only 3x3 in this case, and solve the equations for* using the equations 


x ?r n 



(2p-\) 

X iJ~ 1 


-r M y .* ( 2 v 1) 

1 1 i- 1 ,j 


- U : 


(2p — 2) _ jj 

i+U X i+l,j U iJ + 1 


( 2/7 - 2 ) 

* . 7+1 


(7.25a) 


r ( 2 P) = r _ T (2p-l)_ r (2 P -l) 

x i,/ .y x j-i 


- u 


SW _ 77 


l+ 1 


(2p) 


>2 »+l _ X i,j 


7+1 


(7.25b) 


where p = 1,2, 3, ..., and 


Ki 

= K) 


(7.26a) 

Lij-i 

= b ~j 

L i, 7-1 

(7.26b) 


= B ~ 1 
»>y 

L i~U 

(7.26c) 

V t +ij 

= B~ l 

t J 

U i+ I.j 

(7.26d) 

U iJ+ 1 

= B ~ 1 

* .7 

Ui.j+1 

(7.26e) 


Even though the Jacobian is not updated, in general, for each Newton iteration, the residual vector 
is updated after each Newton iteration. This means that for each Newton iteration a new b must be 
obtained from Eq. (7.26a). This, however, involves only a matrix-vector multiply and not a matrix- 
matrix multiply as is required for L and U, which do not have to be recomputed. By solving Eqs. 
(7.25a) and (7.25b) in place of Eq. (7.23a) and (7.23b) a saving of from zero to a factor of two is 
obtained depending on the number of Newton iterations and symmetric block Gauss-Seidel itera- 
tions that are used. Part of this saving is due to the way the B ’s are inverted. The approach presently 
used is to first perform a decomposition of B using Doolittle’s method. Then invert the two triangu- 
lar matrices resulting from the decomposition, and finally obtain the inverse of B by multiplying the 
two inverted triangular matrices in reverse order [25]. There is no penalty in storage in using the 
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inverse B method given by Eqs. (7.25a) and (7.25b) compared to using Eqs. (7.23a) and (7.23b). 
Experiences gained thus far, indicate that solving Eqs. (7.25a) and (7.25b) is the preferred approach 
for this two-dimensional problem. 
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VIII. Results 


The numerical approach presented here differs somewhat from that presented in [1, 2, and 11]. 
The differences include: (1) the numerical flux formulation at a cell face, (2) the boundary condi- 
tions, (3) as a consequence of (1) the numerical Jacobian is different, and moreover, it is different 
from what one would probably consider to be the customary approach, and (4) the use of the inverse 
of the diagonal block matrices in the solution of the linear system. Numerical results addressing each 
of these changes are given that illustrate the present approach yields results as good or better as be- 
fore, at a reduced computational cost. In addition, a time-dependent solution will be presented and 
compared to a known theoretical solution to demonstrate the accuracy of the method for unsteady 
computations. 

All numerical solutions presented are third-order accurate in space which corresponds to 
x = 1 /3 in the construction of the numerical flux vector as discussed in Chapter V. With the excep- 
tion of the time-dependent solution, all numerical results were carried out for aNACA 0012 airfoil 
with a 161x35 O-type grid. Local time stepping was used for all the steady state solutions. Also, 
all numerical computations were carried out on an IBM RS/6000 Model 560 in 64-bit arithmetic 
using the compilation command “xlf -qautodbl=dblpad -O -bmaxdata:0x3e800000”. 


8.1 Numerical Flux 

The numerical flux formulation used in the past [2] was never given a name. In formulation, 
it falls somewhere between flux extrapolation [26] and MUSCL [9] as explained in [8 and 11]. A 
comparison of the MUSCL approach used here and the no-name flux formulation [8] are compared 
in Fig. 3. For this incompressible flow computation, one will notice that the two flux formulations 
lead to indistinguishable results (except for the region extremely close to the trailing edge where 
MUSCL produce a slightly larger (more positive) pressure). In compressible flow where limiters 
are used, this close agreement between the two numerical flux formulations is not the case [8]. The 
incompressible results given in Fig. 3 are typical of the results obtained to date. 


8.2 Boundary Conditions 

The modification presented in Chapter VI is actually only slightly different from the way bound- 
ary conditions have been treated in the past [1 and 2]. The difference in the two methods is that pre- 
viously phantom points were used and the boundary conditions were determined on the boundaries 
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using characteristic variable boundary conditions and then extrapolated to the phantom points [12]. 
In the new method, the boundary conditions are still determined on the boundaries using the same 
characteristic variable approach, but they are not extrapolated to phantom points, rather they are re- 
tained on the boundaries. This slight modification had an influence on the results in the following 
sense. The converged solutions of the two approaches are essentially indistinguishable as shown 
in Fig. 4, but the convergence rate was changed. Consider the top two curves in Fig. 5 represented 
by the dashed curve and long-dashed curve. The dashed curve is the result from using the old bound- 
ary conditions of [ 1 ] , where ISGS is the number of symmetric block Gauss-Seidel iterations, IFREQ 
represents how often the flux Jacobians are updated, and CFL is, of course, the CFL number. Note 
there is some improvement in the convergence rate brought about by using this slight modification 
in the boundary condition treatment. This new boundary condition treatment however, allows the 
CFL number to be increased considerably as illustrated by the solid curve in Fig. 5. This increase 
comes at the expense of having to increase ISGS, in this case from 5 to 15. This, unfortunately, in- 
creases the CPU time required per iteration. Notice also that the flux Jacobians were updated after 
every iteration (IFREQ= 1 ) with the new boundary conditions for a CFL number of 1000. This also 
increased the CPU time relative to updating them infrequently. However, even with this increase 
in computational work per cycle the end result is a savings in CPU time. For example, the solid line 
in Fig. 5 takes 79 CPU seconds to reach machine zero, whereas the top dashed line in Fig. 5 takes 
109 CPU seconds to reach machine zero. The updating of the flux Jacobians every cycle is, in gener- 
al, not required. By updating the flux Jacobians every cycle for the first 10 cycles and then updating 
every twenty cycles afterwards, the convergence rate of the solid curve with CFL= 1 000 is essentially 
the same and it takes only 38 CPU seconds to reach machine zero. The dash-dot curve in Fig. 5 for 
the old boundary condition treatment with ISGS=15, IFREQ=1, and CFL=40 is included to show 
that the convergence rate is improved relative to using ISGS=5 and IFREQ=20. However, even with 
this extra computational work due to updating Jacobians every cycle, the CFL number could not be 
increased with the old boundary condition treatment as it could with the new boundary condition 
treatment. The CFL numbers used in Fig. 5 are near optimum for each boundary condition treat- 
ment. 

Obviously the new method of handling the boundary conditions is an improvement over the pre- 
vious method, and the improved convergence brought about by the capability to handle a larger CFL 
number is attractive. However, the real pay-off for the capability to handle large CFL numbers is, 
oddly enough, probably greater for unsteady flow than it is for steady state flow. The reason for this 
is that it is always desirable to have the time step restricted by the physics of the problem and not 
the numerics of the solver. This point is discussed in some detail in [27]. 
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Besides increasing the CFL number and improving the convergence rate of the numerical solu- 
tions, this new boundary condition treatment has recently been beneficial in another regard. In using 
a Navier— Stokes version of this numerical approach [1] to compute the flow in the respiratory system 
where there is the problem of handling branching, Gatlin has encountered the problem of the flow 
oscillating considerably among the branches during the transient process of obtaining a steady state 
solution. By using this new boundary condition treatment at the outflow boundary only, the solution 
process was accelerated considerably [28]. 

8.3 Numerical Jacobian 

The solution matrix is usually the Jacobian of a first-order numerical flux vector, even when a 
higher— order numerical flux vector is used in the residual vector. The reason for this is to reduce 
the size of the solution matrix and consequently save storage and usually CPU time. The Jacobian 
used here is developed in Section 7.2. It is not constructed from a first-order numerical flux (unless 
a first-order numerical flux is used in the residual), rather it takes into account some flavor of the 
higher-order numerical flux vector. The reason for trying this approach was due to simplicity and 
reduced operation count, whereas, the reason for keeping this approach is due to the fact it has been 
successful. In fact, the convergence rates experienced in all cases considered thus far have been su- 
perior to the convergence rates obtained with a Jacobian corresponding to a first-order numerical 
flux vector. To illustrate this, consider Fig. 6 which is a comparison of results from the method used 
to obtain the discretized Jacobian presented in Section 7.2 with the results from using the discretized 
Jacobian obtained from a first-order numerical flux vector. Obviously the convergence rate of the 
present method is superior. In this case improved results were obtained at a reduced operation count. 
The reason for this is that a first-order numerical flux does not have to be computed in addition to 
the higher-order numerical flux in order to obtain the discretized Jacobian. 


8.4 Inverse of Diagonal Blocks 

This method of solving the linear system of equations was developed in Section 7.3. In principle 
it should produce the same numerical results as solving the equations without taking the inverse of 
the diagonal blocks as has been the practice in the past [11]. All numerical solutions have been the 
same whether or not the inverse is taken. However, there has been some differences noticed in the 
residual histories when the numbers become small, but this might be due to the different arithmetic 
operations involved in obtaining a solution and the 64-bit accuracy of the computation. 

Timing results of numerical experiments carried out for different values of ISGS (number of 
symmetric block Gauss-Seidel iterations) and IFREQ (number of cycles taken before the numerical 
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Jacobians are updated) are presented in Fig. 7. The times in Fig. 7 are for 100 cycles and all of the 
solutions from which the results of Fig. 7 were extracted were not necessarily converged. These 
results were obtained simply to illustrate the difference in CPU time brought about by using or not 
using the inverse of the diagonal blocks. Three values of IFREQ were considered in Fig. 7. Updating 
the Jacobians every time step corresponds to IFREQ=1 and, of course, this is the most computation- 
ally intense case. For IFREQ= 1 it does not pay to use the inverse of the diagonal blocks for an ISGS 
less than 3 as illustrated by the fact that the solid and dashed curves, with circles, cross at this point 
in Fig. 7. However, for an IFREQ of 9999 (since these solutions were run for 100 cycles an IFREQ 
of 9999 simply means that the Jacobians were computed at the first cycle and then never recomputed 
(frozen) thereafter) Fig. 7 indicates it pays to use the inverse of the diagonal blocks for all values 
of ISGS. Also included in Fig. 7 is a set of curves for IFREQ=5 where, like for IFREQ=9999, it 
pays to use the inverse for all values of ISGS. Figure 7 illustrates that by freezing the Jacobians for 
as few as 5 steps (IFREQ=5) there is a large savings in CPU time compared to updating the Jacobians 
every time (IFREQ=1). Typically steady state solutions are run with an IFREQ of 20. A curve is 
not shown for this, but it is only a small increase above the IFREQ=9999 curves. An IFREQ of 10 
or 20 seems to be sufficient, even for Navier-Stokes computations of complex flows including sepa- 
ration, with extremely tight grids (like maximum y-plus values of 1 and less [29]). For steady state 
solutions of complex problems it is usually not required to update the Jacobians very often after the 
first few cycles as mentioned in Section 8.2 above, so the times corresponding to the IFREQ=9999 
(Jacobian frozen) curves can be achieved in practice. 

Obviously an advantage of taking the inverse of the diagonal blocks is that the number of sym- 
metric block Gauss-Seidel iterations can be increased without as much increase in CPU time as that 
required by not taking the inverse. This, of course, is illustrated in Fig. 7 by the difference in slopes 
of the solid and dashed curves. 

The reader is cautioned that the timings above are for two dimensions, and using the inverse of 
the block diagonals may or may not prove to be as advantageous for three dimensions as it has proven 
to be here for two dimensions. 


8.5 Numerical Solution 

Thus far attention has been focused on speed, storage, and convergence rate. The test case con- 
sidered is rather simple, but the quality of the solution is yet to be demonstrated. To investigate this, 
the numerical solution was compared with a potential flow solution of Bernard [30] in Fig. 8. As 
shown in Fig. 8, the present Euler solution and the potential flow solution [30] are in extremely good 
agreement, except just at the trailing edge of the airfoil where the potential flow solution goes to a 


42 


pressure coefficient of unity whereas the finite volume Euler solution does not. There ,s no cell vol- 
ume centered at the trailing edge of the airfoil in this finite volume formulation and, therefore, the 
finite volume Euler solution does not recompress to a pressure coefficient of exactly unity as the 
node based potential flow solution does. The Euler solution shown in Fig. 8 is considered good. 


8.6 Unsteady Cascade Solution 

To verify an unsteady solution that uses the eigensystem for the time varying curvilinear coordi- 
nate system given in Appendix A, an unsteady solution was obtained for a cascade of infinitely thin 
unstaggered blades oscillating in plunge and in-phase at a reduced frequency based on semi-c or 
of unity This test case was selected because the theory of Smith [31] can be compared to the prese 
numerical solution. The first harmonic is compared with Smith’s theory in Fig. 9, and the compari- 
son is considered good. The grid used was an H-type 121x51 grid with 100 evenly spaced points 
on the blades and the 5 1 points from blade-to-blade were evenly spaced. A complete period of mo- 
tion used 500 time steps and the numerical solution presented was from the fifth period of motion. 
This number of grid points, time steps per period, and number of periods run was over-kill, but the 
objective was to verify the unsteady computational capability of the method and not to try and mini- 
mize any of these parameters. The reason this computation is considered over-kill is based on expe- 
rience in performing computations for this configuration in compressible flow with a similar numer- 
ical approach [32]. 
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IX. Concluding Remarks 


This report describes the most recent computational technology used to solve the two-dimen- 
sional time-dependent incompressible Euler equations. Much of this same technology is used in 
the three-dimensional time-dependent laminar or turbulent, single or multiblock, stationary or dy- 
namic grid, incompressible Navier-Stokes equation codes currently being developed in the CFD 
Lab. Consequently, the material learned will be of use for solving more geometrically and physical- 
ly challenging problems. The report includes considerable detail in certain areas and it is anticipated 
that it might be a helpful aid for teaching some aspects of computational fluid dynamics. 
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Fig. 2 Boundary Condition Terminology 
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COMPARISON OF TWO NUMERICAL FLUX FORMULATIONS 


NACA 0012 At Zero Degrees Angle of Attack 



Fig. 3 Surface Pressure Distributions for the NACA 0012 at Zero Degrees Angle of 
Attack from Two Solutions with Different Numerical Flux Formulations 
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COMPARISON OF TWO BOUNDARY CONDITION TREATMENTS 
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Fig. 4 Surface Pressure Distributions for the NACA 0012 at Zero Degrees Angle of 
Attack from Two Solutions with Different Boundary Condition Treatments 
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Log of the Total Residual 


CONVERGENCE HISTORIES 


NACA 0012 At Zero Degrees Angle of Attack 



Fig. 5 Convergence Histories for Different Boundary Condition Treatments 
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CPU TIME WITH AND WITHOUT BLOCK DIAGONAL INVERSE 


NACA 0012 At Zero Degrees Angle of Attack 



Fig. 7 Computational Times With and Without the Use of the Block Diagonal Inverse 
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Delta Pressure 


CASCADE OSCILLATING IN PLUNGE 


RF=1 , 500 Steps Per Cycle, 121x51 Grid 



Fig. 9 Unsteady Results for a Cascade of Unstaggered Infinitely Thin Blades 
Oscillating in Plunge and In-Phase at a Reduced Frequency of Unity 
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Appendix A 

Eigensystem of the Two-Dimensional Unsteady Incompressible Euler Equations in 
Time-Dependent Curvilinear Coordinates 

The three-dimensional time-dependent incompressible Euler equations in Cartesian coordi- 
nates were transformed to time-dependent curvilinear coordinates in Ref. 1. After the equations 
were transformed to curvilinear coordinates they were written in Chorin’s [3] artificial compressibil- 
ity form as explained in Ref. 1 . Following this same approach, but for the qiore simple case of two 
dimensions rather than three dimensions, the artificial compressibility form of the two-dimensional 
time-dependent incompressible Euler equations in time-dependent curvilinear coordinates is 

dQ , dF_ , dG_ _ q (A. 1) 
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£x = J 1 y v 
= ~ J *1] 
Vx = - J~ l y$ 

Vy = J~ l x i 


and /? = constant. Notice that flux vectors F and G can be written 


where 


K = J 


ft t ^r) 


u6 k + k x p 

vd k + ky p 


0 k = n + v + fc, 


and 


K = F and 6 k = U for k = £ 
K = G and 6 k = V for k = rj 


The quasilinear form of Eq. (A.l) is 


^ + a§ + sf 

dr d£ dt] 


= 0 


where 


A = M. B = £G 

dQ dQ 


(A. 2) 


(A.3) 


(A. 4) 


(A.5) 


The Jacobian matrices A and B are denoted by K where the elements of K are given by 
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(A.6) 
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Therefore, 


and 


K = F and K = A when k = £ 


K — G and K = B when k = rj 


The flux vectors represented by Eq. (A. 2 ) can be written in terms of the elements of the dependent 

T 

variable vector Q = [ Jp , Ju, Jv] T = | Q\ Q 2 , ^3] as 
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(A. 7 ) 


The elements of the Jacobian matrix K are determined from Eq. (A.6) resulting in 
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(A. 8) 
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The eigenvalues of K are given by 
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and 6 k is given by Eq. (A. 3). Therefore, 
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(A. 11) 


where the subscript k corresponds to k = £ or k = rj. The corresponding right eigenvectors of 
K are given by the columns of 
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The corresponding left eigenvectors of K are given by the rows of 
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(A. 13) 


The eigensystem of K is therefore composed of the eigenvalues given by Eqs.(A.9), the right 
eigenvectors given as the columns of the matrix R k (Eq.(A.12)), and the left eigenvectors given as 
the rows of the matrix R k 1 (Eq.(A. 13)). Note that each eigenvector corresponds to one particular 
eigenvalue. The Jacobian matrices are now diagonalized by 

R~ l KR k = A k (A. 14) 
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